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TRIANGULAR  DECOMPOSITION  Of  A  POSITIVE  DEFINITE  MATRIX 
PLUS  A  SYMMETRIC  DYAD  WITH  APPLICATIONS  TO  KALMAN  FILTERING 


i 


I.  INTRODUCTION .  Given  a  positive  definite  matrix  P  Choleski's  theorem 
states  that  there  exists  a  real,  non-singular,  lower  triangular  matrix 
L  such  that 


LLT  =  P  (1)  j 

I 

i 

Furthermore  if  the  diagonal  elements  of  L  are  taken  to  be  positive  the  I 

decomposition  is  unique.  L  is  called  the  square  root  of  F.  The  Choleski 

algorithm  for  decomposition  of  P  is  presented  in  [1]  and  [2].  The  ! 

Choleski  decomposition  is  very  useful  in  many  numerical  linear  algebra 

i 

problems.  In  particular  it  provides  a  useful  numerical  technique  in  the  \ 

matrix  square  root  formulation  of  the  Kalman  filter  [3].  The  triangular  ' 

i 

decomposition  of  Choleski  is  extended  below  to  the  decomposition  of  a  posi-  ! 

T  ! 

tive  definite  matrix  P  plus  a  symmetric  dyad  cxx  . 

i 

II.  FIRST  TRIANGULAR  DECOMPOSITION  ALGORITHM.  Suppose  we  have  a  lower  j 

trianguleir  decomposition  L  of  a  positive  definite  matrix  P.  The  elements 
of  L  must  satisfy  the  equations. 

9  i 

i 

i  ; 

l  *(i,j)Mk,j)  =  P.k  k>i  (2) 

j=i  i 

i 

l  i2(i,j)  =  Pu  (3) 

j=i  i 

i 

\ 

i 


.j 

j 


Consider  the  problem  of  computing  the  triangular  decomposition  b'  of  the 
modified  matrix 


P  =  P  +  cxx 


(4) 


given  the  decomposition  L  of  P.  The  decomposition  b'  must  satisfy  the 
equations 


i  i 

l  X/(i,j)5/(k,j)  =  l  Hi, j)H(k,j)  +  k>i 

j=l  j=l 


(5) 


I  l U,j) 
j=l 


l  «.2(i,j)  +  cx^ 
j=l 


(6) 


S'- 

>.  ■ 


First  consider  (5)  and  (6)  for  the  case  of  i=l.  In  this  case 


K.'(.l,l)i'(k,l)  =  Ul,lU(k,l)  +  cxxxk  k>l 


(7) 


r2(i,l)  =  a2(lfl)  +  cx2 


(8) 


The  first  column  of  the  modified  matrix  b'  is  easily  computed  from  (7)  and 
(8).  Now  rewrite  (C)  as 


i  l 

l  r(i,j)r(k,j)  +  r(i,in'(k,i)  =  i  ui,jmk,j)  +  la.nMk.u 

j=2  j=2 

(9) 


+  cx.x, 
i  k 


2 


The  second  term  on  the  left  of  (9)  can  be  computed  from  (7)  as 


ni.DHk.l) 


1  !^1>1)  A(i,l)Jt(k,J  )  +  ex  x  + 

£'2(1,1)  t'2(l,l) 


(10) 


l(ktl)l(l^l)  +  11 L-  X.x. 
r2u,i)  1  1  It  (1,D  * 


Substituting  (10)  into  (9)  and  combining  terras  gives 


i  i 

l  ru,j)r(k,j)  =  I  Jt(i,j)Jt(k,j)  *  c(iyi)41) 

j=2  J-2 

where 


(1)  _  cft2(ltl) 
c  O 

1  (1,1) 


x1Ui,D\ 

xnrn  J 


Similarly,  (6)  can  be  written  as 


(11) 


(12) 


(13) 


i  i 

X  Jt'2(i,j)  +  i'2(i,l)  =  l  Jl2(i,j)  +  A2(i,l)  +  cx2  (14) 

j=2  j=2 
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Substituting  from  (7)  for  the  second  ter™  on  the  loft  of  (14)  And  combining 
terms  gives 

i  i  2 

l  *'2U,1>  =  I  *2U,j)  +  c(1)x<ii)  (1&) 

j=2  j=2 

Equations  (11)  and  (15)  define  a  decomposition  problem  equivalent  to  the 

original  problem  defined  by  (5)  and  (6)  but  with  the  dimension  reduced 

by  one.  It  is  easily  seen  that  the  application  of  the  alxovo  technique 

n  times,  each  time  reducing  the  dimension  of  the  decompop it  ion  problem  by 

one,  solves  the  original  decomposition  problem.  The  following  equation* 

T 

summarize  the  algorithm  for  tnc  triangular  decomposition  of  f’+cxx  . 


III.  SECOND  TRIANGULAR  DECOMPOSITION  ALGORITHM.  An  alternate  decomposi¬ 
tion  of  positive  definite  matrix  is  possible.  Let  L  be  i  unit  lower  tri¬ 
angular  matrix,  i.e.,  having  ones  along  the  diagonal.  Thun  a  poiiitive 
definite  matrix  P  can  be  written  as 


4 


P  =  LDLT  •  (20) 

where  D  is  a  diagonal  matrix  of  positive  numbers.  An  algorithm  for 
computing  the  decomposition 

,  0. 

L'd'L'T  =  LDLT  +  cxxT  (21) 

can  be  derived  in  a  manner  parallel  to  the  previous  decomposition 
algorithm.  The  algorithm  is  as  follows 


d^ 

l 


♦  c(i>x<«‘ 
1 


i  =  1  ,n 


(22) 


^x+D  =  x(i)  _  x(1^(k#i) 


(23) 


£^(k,i)  ~  &(k,i)  + 


c(i)x<» 


(  i+1 ) 
*k 


} 


i  +  l<k<n 


c(itl) 


(24) 


(25) 


IV,  APPLICATION  TO  KALMAN  FILTERING  I.  At  a  measurement  update  in  a 
discrete  Kalman  filter  the  state  estimate  io  given  by 


x 


k/k-I 


pAl*ki,Vi) 


(26) 


5 


where  the  covariance  matrix  P.  is 

k 

pk  =  (p <”> 

Consider  the  measurement  z,  to  be  a  scalar  since  the  vector  measurement 

k 

problem  can  always  be  reduced  to  this  case,  see  [33.  In  this  scalar  case 
is  a  row  vector,  is  a  scalar,  and  is  n*n.  In  the  matrix  square 
root  formulation  of  the  Kalman  filter  given  in  [33,  (27)  is  the  inversion 
of  a  matrix  plus  a  dyad  which  is  simply  expressed  by  the  method  of  modi¬ 
fication  given  in  Householder  [4].  This  gives 


Pk/k-l 


?k/k-lW*k/k-l 


(28) 


since 


Pfc  and  ?  .^  are  positive  definite  write  them  as 


and 


k/k-1 


Lk/k-lLk/k-l 


Anoivne  tlujt  Lyjy  ^  ie  lower  triangular  with  positive  diagonal  elements. 
Gubol ituting  in  (28) 


T 

L. 

k  k 


Jk/k-J 


T  T  T 

v  -1  ’  ckLk/k-lLk/k-lHkHkLk/k-lLk/k-l 


(29) 


=  3/(  vvwX ) 


(30) 


Let 


“k  = 


and 


Wk  =  Lk/k-lUk 


(31) 


(32) 


Then 


Vv  *  Lk/k-iLk/k-i  -  wl  (33 

Thus  (33)  is  in  the  proper  form  for  application  of  the  first  algorithm. 
In  addition  to  the  above  computations  the  state  estimate  given  by  (26) 
must  also  be  computed.  Some  manipulation  shows  that  (26)  can  be  written 
as 


*k  = 


k/k-1 


+  c 


cWk(  VHkxk/k-l ) 


(34) 


Now  consider  the  numerical  efficiency  of  the  application.  The 
execution  of  (31)  and  (32)  takes  n(n+X)  (mult).  (30)  takes  n  (mult)  and 
1  (div),  (34)  takes  (n+1 )  (mult).  The  decomposition  algorithm  for  03) 
requires 

(3n2+9n-6)  /2  (mult) 
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2(n-i)  (div)  and  n  v^s.  Thus  the  use  of  the  first  decomposition  algo¬ 
rithm  requires 

^ S/2  n2+15/2  n-2 j  (mult) 

2n-l  (div)  and  n  /"'s.  Tha  technique  presented  in  [33  which  we  have  been 

2 

using  in  our  Kalman  filter  program  requires  about  3n  +2n  (mult)  but  does 
not  generate  a  triangular  square  root  matrix. 

Now  consider  the  numerical  efficiency  of  the  second  algorithm.  A 
development  paralleling  (29)-(34)  gives 


where 


LkDkI*k  =  Sc/k-l^/k-lh^/k-l  ‘  ckWkWk 

(35) 

Uk=  Vk-lS/k-A 

(36) 

wk  =  Lk/k-l\ 

(37) 

‘k  * i-  ( Wk/k-X ) 

(38) 

The  use  of  the  second  algorithm  and  execution  of  (34)-(38)  requires 
2 

2n  +7n  (mult)  and  n  (div). 


The  application  of  the  second  algorithm  results  in  fewer  operations 
than  either  our  present  algorithm  or  the  first  algorithm  presented 
above  plus  the  benefit  of  having  a  covariance  square  root  which  is 
triangular. 


V*  APPLICATION'  TO  KALMAN  FILTERING  II.  Rather  than  compute  the  square 
root  of  the  covariance  matrix  the  square  root  of  the  inverse  covariance 
matrix  may  be  computed.  The  updating  of  the  inverse  covariance  matrix 
is  given  by  (27). 


-  P 


-1 

k/k-1 


+ 


VkX 


(39) 


since  the  inverse  covariance  is  also  positive  definite,  write 


p'1  =  l-1  d_1  rT 

k/k-1  k/k-1  uk/k-l  Lk/k-l 


and 


skA  ■ 


T 

Again  considering  scalar  measurements  let  Uj,  =  and  ck  =  1/R^.  Then 


k/k-lDk/k-lLk/k-l  +  ckukuk 


(40) 


which  is  in  the  form  for  application  of  the  second  algorithm.  The 
computation  of  the  corrected  state  estimate  given  by  (26)  requires  the 
computation  of 


y«,  = 


P_X 


(  zk'Hkxk/k-l  ) 


(41) 


I 

I 


i 

t 

i 

t 

> 

i 


| 


s 

i 

: 


I 

i 


i 
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Let  e,  =  c. 
k  k 


(  VHk*k/k-l  ) 


and  rewrite  (41)  as 


Pklyk  =  Lk\LkTyk  =  9k 


The  solution  of  (42)  for  y^  requires  the  solution  of  two  sets  of  linear 

equations  each  with  a  triangular  coefficient  matrix.  The  solution  of  a 

triangular  set  cf  linear  equations  is  a  standard  procedure  in  numerical 

2 

analysis.  The  solution  of  (42)  for  the  vector  requires  about  n  +n 
(mult)  and  n  (div).  The  solution  of  the  Kalman  filter  equations  (40) 
and  (42)  at  a  measurement  update  requires  about  n  +5n  (mult)  and  n  (div) 
using  the  second  decomposition  algorithm.  The  specification  of  (40)  and 
(42)  as  the  measurement  update  equations  requires  that  l  comPut®<l 

at  a  time  update.  In  the  special  case  where  there  are  no  process 
dynamics,  i.e.  we  are  estimating  constant  parameters,  and  there  is  no 
process  noise  Lj^.j/L”^  50  that  no  additional  computation  is  required 
to  obtain  L~^  Although  this  is  a  special  case,  it  is  of  interest  in  the 
bias  filter  portion  of  the  WSHR  BET  program  where  the  measurement  biases 
are  assumed  to  be  constant  in  time. 


VI.  APPLICATION  TO  KALMAN  FILTERING  III.  The  use  of  the  second  tri¬ 
angular  decomposition  algorithm  in  the  WSMR  BET,  which  is  an  extended 
Kalman  filter,  has  resulted  in  a  significant  increase. in  numerical 
efficiency,  however,  the  motivating  factor  for  the  development  and  use 
of  this  algorithm  was  for  the  application  described  below. 

The  Kalman  filter  in  the  WSMR  BET  is  divided  into  two  filters,  the 

* 

zero-bias  filter  which  produces  trajectory  state  estimates  x^  and  the 


10 


bias  filter  which  produces  estimates  of  the  measurement  biases.  These 
two  estimates  are  combined  to  form  the  optimal  trajectory  state  estimate 
defined  by 


*k  = 


T  b 
k  k 


(43) 


where  is  a  combining  matrix  which  must  satisfy  equations  determined  from 

the  orthogonality  properties  of  the  statr-  estimates.  One  property  of  the 

.  *  A 
filter  decomposition  given  in  (43)  is  the  requirement  that  and  b^  be 

orthogonal.  This  requirement  is  naturally  satisfied  by  the  equations 
governing  x^,  b^,  and  T^  except  at  a  point  where  a  measuring  instrument 
is  deleted  from  the  filtering  process.  An  instrument  may  be  dropped  be¬ 
cause  it  is  no  longer  taking  observations,  its  bias  is  too  large,  or  the 
measurements  are  chronically  inconsistent  with  their  statistics.  Let 

*  A 

x  ,  P  denote  the  zero-bias  state  estimate  and  its  covariance  just  prior 
~  *  * 
to  dropping  a  measurement  from  the  filtering  solution  and  let  x  and  P 

be  the  same  quantities  immediately  after  dropping  the  measurement. 

A  A 

Similarly  let  b  and  bx  denote  the  bias  state  estimates  before  and  after 

—  *r  /\  a 

dropping  a  measurement.  b+  is  formed  by  deleting  the  component  of  b 
corresponding  to  the  measurement  being  dropped.  T  and  T+  are  the 
combining  matrices  before  and  after.  T+  has  one  less  column  than  T  . 

Pfc  and  P^h  are  the  bias  covariance  matrices  before  and  after  dropping  a 
measurement.  P^+  is  formed  from  P^_  by  deleting  the  bow  and  column 
corresponding  to  the  measurement  being  dropped.  The  updating  equation 
for  x  is 


(44) 
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where  t.  is  the  column  of  T  being  deleted,  b,  is  the  bias  estimate  for 

*  *  X  ^ 

the  measurement  being  dropped,  and  l  is  vector  chosen  so  that  x+  and  b+ 
will  be  orthogonal  in  the  usual  statistical  sense.  The  updating  equation 
for  P  is 


ft 

P 

+ 


*' 

P  + 


-M 


b- 


(i,i)-S,TP. 


b+ 


(45) 


If 


ft  ft  ft  ftT 

?  =  C  D  C 


and 


ft 

P 

+ 


ft  ft  *T 
C+DtCt 


then  (45)  is  in  a  form  for  application  of  the  second  triangular  decompo¬ 
sition  algorithm.  In  addition,  let 


cb-Vc 


T 

b-1 


and 


Pb+  "  Cb+Db+Cb+ 


Let  the  i**1  row  and  column  C.  be  deleted.  Also  let  the  i**1  row  and 

b- 

cclumn  of  D.  be  deleted.  Call  the  resulting  matrices  C'  and  D '  . 

b-  o-  b- 


12 


Then 


VW  *  cb-“b-ct 


Diei°i 


(46) 


where  is  the  column  deleted  from  C^_  and  d.  is  the  diagonal  element 
deleted  from  D^_.  Thus  (46)  presents  another  application  for  the 
decomposition  algorithm. 
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